# This file generates figures and tables for the dissection part of the paper
library(HeadR)
library(latex2exp)
setwd("~/Dropbox/BLP_REStat/Replication")
#
rm(list=ls())
#
# Mixed logit ====
mix <- "MLOG/MC"
OG <- 90
for(pan in c("a","b")) {
# content of table  
monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,".rds"))
CG <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),change.ces=100*mean(change.ces),
                           sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                        pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif),eta=mean(eta)
                        ),by=.(v)]
CG[, v_name:="Mixed Logit"]
CG[v==1, v_name:="Logit"]
CG[v==2, v_name:=" $\\beta $ heterogeneity"]
CG[,output := texout(.(v_name,change.mixed,change.ces,eta,pt,pte))]
writeLines(CG$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,".tex"))
#
d <- 0.15 # separating the two bars
#
# changes in domestic market share graph
pdf(file=paste0("Graphs/CF_MM/",mix,"/changes_",pan,"OG_",OG,".pdf"),height=3.5,width=4)
par(mar=c(4,5,1,1))  #bottom, left, top and right
ymax <- CG[,max(change.mixed+2*sd.change.mixed)]
CG[,plot(v+d,change.mixed,type="h",
              ylab="Increase in domestic share \n (percentage points)",xlab="Parameter setting",
              lwd=10,lend=1,ylim=c(0,10),yaxs="i",xlim=c(0.5,3.5),col="grey70",xaxt="n")]
CG[,axis(side=1,at=v,labels=c("Logit",TeX(r'($\beta$ het)'),"Mixed Logit"))]
CG[,lines(v-d,change.ces,type="h",col="grey30",lwd=10,lend=1)]
CG[,arrows(v-d,change.ces+1.96*sd.change.ces,v-d,change.ces-1.96*sd.change.ces,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
CG[,arrows(v+d,change.mixed+1.96*sd.change.mixed,v+d,change.mixed-1.96*sd.change.mixed,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
legend("topleft",col=c("grey30","grey70","black"),lwd=c(10,10,3),bty="n",
       legend=c("Approx Mean","True Mean","95% CI"))
dev.off()
}

# Now graph the setting 3 ("BLP") for OG 90, 50, 10
pan <- "b"
for(OG in c(90,50,10)) {
  monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,".rds"))
  assign(paste0("changes",OG),monte.reps[v==3,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),change.ces=100*mean(change.ces),
                           sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                           pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif),eta=mean(eta)
  ,og=OG),by=.(v)])
}
CG <- rbind(changes90,changes50,changes10)
CG[,v := 1:3]
pdf(file=paste0("Graphs/CF_MM/",mix,"/changes_",pan,"_allOG",".pdf"),height=3.5,width=4)
par(mar=c(4,5,1,1))  #bottom, left, top and right
ymax <- CG[,max(change.mixed+2*sd.change.mixed)]
CG[,plot(v+d,change.mixed,type="h",
              ylab="Increase in domestic share \n (percentage points)",xlab="Outside good share (%)",
              lwd=10,lend=1,ylim=c(0,12),yaxs="i",xlim=c(0.5,3.5),col="grey70",xaxt="n")]
CG[,axis(side=1,at=v,labels=og)]
CG[,lines(v-d,change.ces,type="h",col="grey30",lwd=10,lend=1)]
CG[,arrows(v-d,change.ces+1.96*sd.change.ces,v-d,change.ces-1.96*sd.change.ces,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
CG[,arrows(v+d,change.mixed+1.96*sd.change.mixed,v+d,change.mixed-1.96*sd.change.mixed,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
legend("topleft",col=c("grey30","grey70","black"),lwd=c(10,10,3),bty="n",
       legend=c("Approx Mean","True Mean","95% CI"))
dev.off()
#
CG[,output := texout(.(og,change.mixed,change.ces,eta,pt,pte))]
writeLines(CG$output,paste0("Tables/CF_MM/",mix,"/changes2.tex"))

# Mixed CES ====
mix <- "MCES/MC"
OG <- 90
pan <- "b"
monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,".rds"))
changes <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),change.ces=100*mean(change.ces),
                         sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                         pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif),eta=mean(eta)
),by=.(v)]
changes[, v_name:="Mixed CES"]
changes[v==1, v_name:="CES"]
changes[v==2, v_name:=" $\\beta $ heterogeneity"]

changes[,output := texout(.(v_name,change.mixed,change.ces,eta,pt,pte))]
writeLines(changes$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,".tex"))
#
d <- 0.15 # separating the two bars
#
# Changes in domestic market share
pdf(file=paste0("Graphs/CF_MM/",mix,"/changes_",pan,"OG_",OG,".pdf"),height=3.5,width=4)
par(mar=c(4,5,1,1))  #bottom, left, top and right
ymax <- changes[,max(change.mixed+2*sd.change.mixed)]
changes[,plot(v+d,change.mixed,type="h",
              ylab="Increase in domestic share \n (percentage points)",xlab="Parameter setting",
              lwd=10,lend=1,ylim=c(0,10),yaxs="i",xlim=c(0.5,3.5),col="grey70",xaxt="n")]
changes[,axis(side=1,at=v,labels=c("CES",TeX(r'($\beta$ het)'),"Mixed CES"))]
changes[,lines(v-d,change.ces,type="h",col="grey30",lwd=10,lend=1)]
changes[,arrows(v-d,change.ces+1.96*sd.change.ces,v-d,change.ces-1.96*sd.change.ces,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
changes[,arrows(v+d,change.mixed+1.96*sd.change.mixed,v+d,change.mixed-1.96*sd.change.mixed,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
legend("top",col=c("grey30","grey70","black"),lwd=c(10,10,3),bty="n",
       legend=c("Approx Mean","True Mean","95% CI"))
dev.off()

#
mix <- "MCES/MC"
OG <- 10
pan <- "b"
monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,".rds"))
changes <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),change.ces=100*mean(change.ces),
                         sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                         pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif),eta=mean(eta),
                         avg.diff=100*mean((change.ces-change)),se.diff=100*sd((change.ces-change))/sqrt(.N)),by=.(v)]
changes[, v_name:="Mixed CES"]
changes[v==1, v_name:="CES"]
changes[v==2, v_name:=" $\\beta $ heterogeneity"]
changes[,output := texout(.(v_name,change.mixed,change.ces,eta,pt,pte))]
writeLines(changes$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,".tex"))
changes.b <- copy(changes)

#
mix <- "MCES/OLY"
OG <- 10
pan <- "c"
monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,".rds"))
changes <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),change.ces=100*mean(change.ces),
                         sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                         pt=mean(pt),pte=mean(pte),reg.shr.dif=round(mean(reg.shr.dif),3),eta=round(mean(eta),3),
                         avg.diff=100*mean((change.ces-change)),se.diff=100*sd((change.ces-change))/sqrt(.N)),by=.(v)]
changes[, v_name:="Mixed CES"]
changes[v==1, v_name:="CES"]
changes[v==2, v_name:=" $\\beta $ heterogeneity"]
changes[,eta := eta +0.00001] # add a small number within tolerance
changes[,output := texout(.(v_name,change.mixed,change.ces,eta,pt,pte))]
writeLines(changes$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,".tex"))
changes.c <- copy(changes)

# Graphs of predicted and true changes in domestic market share
d <- -0.15 # separating the two bars
pdf(file=paste0("Graphs/CF_MM/",mix,"/diffs_OG_",OG,".pdf"),height=3.5,width=4)
par(mar=c(4,5,1,1))  #bottom, left, top and right
yrng <- as.vector(changes.b[,.(min(avg.diff-2*se.diff),max(avg.diff+2*se.diff))])
changes.b[,plot(v+d,avg.diff,type="h",
              ylab="Bias in domestic share changes \n (percentage points)",xlab="Parameter setting",
              lwd=10,lend=1,ylim=c(-1.2,0.8),yaxs="i",xlim=c(0.5,3.5),col="grey30",xaxt="n")]
changes.b[,axis(side=1,at=v,labels=c("CES",TeX(r'($\beta$ het)'),"Mixed CES"))]
abline(h=0)
changes.c[,lines(v-d,avg.diff,type="h",col="grey70",lwd=10,lend=1)]
changes.c[,arrows(v-d,avg.diff+1.96*se.diff,v-d,avg.diff-1.96*se.diff,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
changes.b[,arrows(v+d,avg.diff+1.96*se.diff,v+d,avg.diff-1.96*se.diff,col="black",
                lwd=3,angle=90,length=0.02,code=3)]
points(1-d,0,pch=15,col="grey70")
legend("bottomleft",col=c("grey30","grey70","black"),lwd=c(10,10,3),bty="n",
       legend=c("Monop. Comp.","Oligopoly EHA","95% CI"))
dev.off()



# New tables including results using eta estimated on aggregate data

# Mixed CES ====
mix <- "MCES/MC"
OG <- 90
pan <- "b"
monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,"_aggeta.rds"))
CG <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),
                    change.ces=100*mean(change.ces),change.ces.agg=100*mean(change.ces.agg),
                    sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                    pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif)
                    ,eta=mean(eta),aggeta=mean(aggeta)
),by=.(v)]
CG[, v_name:="Mixed CES"]
CG[v==1, v_name:="CES"]
CG[v==2, v_name:=" $\\beta $ heterogeneity"]
#
CG[,output := texout(.(v_name,change.mixed,change.ces,change.ces.agg,eta,aggeta,pt,pte))]
writeLines(CG$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,"_aggeta.tex"))



# Mixed logit ====
mix <- "MLOG/MC"
OG <- 90
pan <- "b"
  monte.reps <- readRDS(file=paste0("Tables/CF_MM/",mix,"/monte_allv_",pan,"OG_",OG,"_aggeta.rds"))
  CG <- monte.reps[,.(dom.shr=100*mean(dom.shr),change.mixed=100*mean(change),
                      change.ces=100*mean(change.ces),change.ces.agg=100*mean(change.ces.agg),
                      sd.change.mixed=100*sd(change), sd.change.ces=100*sd(change.ces),
                      pt=mean(pt),pte=mean(pte),reg.shr.dif=mean(reg.shr.dif)
                      ,eta=mean(eta),aggeta=mean(aggeta)
  ),by=.(v)]
  CG[, v_name:="Mixed Logit"]
  CG[v==1, v_name:="Logit"]
  CG[v==2, v_name:=" $\\beta $ heterogeneity"]
  
  CG[,output := texout(.(v_name,change.mixed,change.ces,change.ces.agg,eta,aggeta,pt,pte))]
  writeLines(CG$output,paste0("Tables/CF_MM/",mix,"/changes1",pan,"OG_",OG,"_aggeta.tex"))
  

